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One goal of computational anatomy (CA) is to develop tools to accurately segment 
brain structures in healthy and diseased subjects. In this paper, we examine the 
performance and complexity of such segmentation in the framework of the large 
deformation diffeomorphic metric mapping (LDDMM) registration method with reference 
to atlases and parameters. First we report the application of a multi-atlas segmentation 
approach to define basal ganglia structures in healthy and diseased kids' brains. The 
segmentation accuracy of the multi-atlas approach is compared with the single atlas 
LDDMM implementation and two state-of-the-art segmentation algorithms — Freesurfer 
and FSL— by computing the overlap errors between automatic and manual segmentations 
of the six basal ganglia nuclei in healthy subjects as well as subjects with diseases 
including ADHD and Autism. The high accuracy of multi-atlas segmentation is obtained at 
the cost of increasing the computational complexity because of the calculations necessary 
between the atlases and a subject. Second, we examine the effect of parameters on total 
LDDMM computation time and segmentation accuracy for basal ganglia structures. Single 
atlas LDDMM method is used to automatically segment the structures in a population of 
16 subjects using different sets of parameters. The results show that a cascade approach 
and using fewer time steps can reduce computational complexity as much as five times 
while maintaining reliable segmentations. 
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INTRODUCTION 

Computational anatomy (CA) methods have been at the forefront 
of neuroimaging studies of neurodevelopment and neurodegen- 
eration. Software such as FSL (Jenkinson et al., 2012), AFNI 
(Cox, 2012), BrainVoyager (Goebel, 2012), Caret (Van Essen, 
2012), FIASCO (Aguirre, 2012), SUMA (Saad and Reynolds, 
2012), FreeSurfer (Fischl, 2012), and SPM (Ashburner, 2012) 
have been used to analyze structural and functional properties of 
the human brain via magnetic resonance images (MRI), diffusion 
tensor images (DTI), and functional MRI (fMRI) at a resolution 
of 1 mm 3 . In addition, pipelines for registering and visualiz- 
ing hundreds of images have emerged. These include DTIStudio 
(Jiang et al, 2006), LONI (Dinov et al., 2006), ANTs (Avants 
et al., 2008), and 3D Sheer (Pieper et al., 2004). At the heart 
of these methods is the modeling of anatomical structures and 
their huge variation, an active field originally formalized as CA 
(Grenander and Miller, 1998). The main difficulty is the com- 
plexity of anatomical structures and the large variation between 
individuals. Here anatomical structures are represented as a col- 
lection of coordinate systems: landmark points, curves, surfaces, 



and sub-volumes. These structures are represented as deformable 
templates, with the space of anatomical images being the set gen- 
erated by the group of diffeomorphic transformations acting on 
the template with associated probability laws, which describe their 
variation. The transformations are detailed so that a large family 
of shapes may be generated with the precise topology of the tem- 
plate maintained. In particular, our methods have demonstrated 
localized shape differences in multiple sub-cortical structures in 
neuroimaging studies of Alzheimer's Disease (Qiu et al., 2008a, 
2009c), ADHD (Qiu et al., 2009b), Autism (Qiu et al., 2010), 
schizophrenia (Wang et al., 2008), and Tourette Syndrome (Wang 
et al., 2007). 

The theoretical framework we have adopted is based on the 
large deformation diffeomorphic metric mapping (LDDMM) 
algorithm (Beg et al., 2005) and advances have been developed 
by others as well (Risser et al., 2010, 2011; Auzias et al., 2011). 
It also has been demonstrated that the best registration methods 
have high dimensional and diffeomorphic properties that have 
in common many of the properties incorporated into LDDMM 
(Klein et al., 2009, 2010). The LDDMM approach allows shape 
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to be uniquely encoded by the vectors normal to the outline 
of the template (Miller et al, 2006). This property provides the 
basis of our shape analysis projects and is thus crucial in locat- 
ing changes in multiple structures caused by disease leading to 
a better understanding of the effect of disease on neighboring 
structures. 

To date, many MRI studies of subcortical gray matter nuclei 
have defined a single measure of structural volume. While this has 
the advantage of being quantitative, it does not give specific infor- 
mation about subregions of atrophying nuclei. This information 
would be useful in order to determine whether MRI morphomet- 
ric results correlate with neuropathologic studies, to define better 
the subregional distribution of atrophy, and for correlation of 
pathologic changes with clinical features of the disease. Methods 
of statistical shape analysis have proved useful for studying nor- 
mal age related changes in subcortical nuclei and for studying a 
number of other diseases (Bansal et al., 2007; Brun et al., 2009; 
Chung et al., 2010; Stein et al., 2011; van den Bogaard et al., 
2011). Such analysis based on our LDDMM framework (Miller 
and Qiu, 2009) has proved useful for studying changes in sub- 
cortical nuclei in normal aging and several diseases (Csernansky 
et al, 1998, 2000, 2002; Wang et al, 2007; Qiu and Miller, 2008; 
Qiu et al, 2008b, 2009a,b,c,d, 2010). For example, analysis with 
the LDDMM based surface based morphometry (SBM) pipeline 
revealed heterogeneity of atrophy in subcortical structures in 
studies of Huntington's disease (Younes et al., 2012) and dementia 
(Miller etal., 2012). 

So as the usage of neuroimaging software expands, it becomes 
critically important to provide novel analyses that go beyond 
volumes and/or other scalar quantities at the voxel or regional 
level. Specifically, the ability to segment subcortical structures in 
hundreds of scans which can then be used for shape analysis. 
However, segmentation accuracy and performance is limited by 
two complementary factors. 

The first is the restriction to a single atlas as in MRIStudio 
which has become widely used in neuroimaging studies of MRI 
and DTI data as evidenced by meta-data search in Google Scholar. 
Currently, MRIStudio consists of three programs: DTIStudio 
(Jiang et al, 2006), DiffeoMap and ROIEditor. In particular, 
DiffeoMap is a program for image transformation based on 
LDDMM and ROIEditor uses the results of DiffeoMap to per- 
form image analysis with respect to a single atlas both at the 
voxel and regional level. The restriction can be compounded by 
the wide anatomical variability in structures. For example, ven- 
tricles differ in shape and geometry in neuroimaging studies of 
aging and dementia (Reig et al., 2007). Further, proximity of 
sub-cortical structures can affect the accuracy and performance 
of automated segmentation algorithms resulting in conflicting 
results at the volumetric level (Heckemann et al., 2010). Here we 
explore the incorporation of multi-atlases to drive the segmen- 
tation of subcortical structures which should lead to improved 
statistical significance in shape analysis even from a small sample. 

The segmentation problem, which relies on learning anatom- 
ical information from pre-segmented training datasets, is usually 
handled in the Bayesian framework by solving a maximum a pos- 
teriori (MAP) estimation problem. The information of multiple 
pre-segmented training datasets can be combined naturally in the 



Bayesian framework. Compared with single atlas based segmenta- 
tion, multi-atlas based segmentation has been shown to be more 
powerful and more accurate (Rohlfing et al., 2004; Heckemann 
et al, 2010; Langerak et al, 2010; Lotjonen et al, 2010). More 
recently we have demonstrated the use of LDDMM with multi- 
ple atlases in which both the diffeomorphic coordinate change as 
well as the atlas being selected is unknown. Such an approach has 
been shown to be efficient in mediating large deformations in the 
context of building reliable segmentations of both subcortical and 
cortical structures (Tang et al., 2012). 

This leads to the second problem which is the computational 
complexity of LDDMM which can be compounded by multi- 
ple atlases. By default, LDDMM computes the geodesic along 
the Riemannian manifold from the template to the target in 
ten uniform time steps (Beg et al., 2005; Ceritoglu et al., 2009; 
Tward et al., 2011). As such it takes about 30min to map one 
whole brain at 1 mm 3 resolution using a parallelized implemen- 
tation on Intel Xeon CPU E5530 at 2.40GHz with 16 cores. By 
evolving the flow field over time, it allows for large deforma- 
tions which come at a cost in terms of computational resources. 
In the multiple atlas setting, the computational requirements 
increase as order number of atlases. In an attempt to overcome the 
computational complexity, LDDMM has inspired alternative dif- 
feomorphic implementations image registration algorithms such 
as DARTEL (Ashburner, 2007), Diffeo DEMONS (Vercauteren 
et al., 2009) and Spherical Demons (Yeo et al., 2010) to name 
but a few. These methods essentially compute the mapping in 
one time step to generate a "stationary velocity field" in the spirit 
of exponential maps for finite dimensional groups. While these 
mappings no longer satisfy the conservation law and therefore 
metric property giving rise to complete encoding of the diffeo- 
morphic flow of the shape in the initial tangent vector (Miller 
et al., 2006), their diffeomorphic features make them power- 
ful alternatives for generating segmentations which are accurate. 
Then metric structure can be derived from the segmentations as 
demonstrated (Miller and Qiu, 2009). By reformulating the opti- 
mization problem, LDDMM can be recast as EPDiff, i.e., an initial 
value problem (IVP) in which geodesic shooting is used to evolve 
the initial momenta to match with the target shape (Miller et al., 
2006; Marsland and McLachlan, 2007; Ashburner and Friston, 
2011; Vialard etal, 2012). 

This paper examines performance and complexity of LDDMM 
with reference to atlases and parameters. A natural way is to 
examine tradeoff of computational complexity vs. segmentation 
reliability. Hence we examine segmentations of multiple sub- 
cortical structures from ongoing large scale neuroimaging studies 
which are different and more complex than those previously 
studied (Tang et al, 2012, 2013). 

METHODS AND DATA 
DATA 

The first dataset included thirty whole brain, high resolution Tl- 
weighted 3D-volume MPRAGE images (matrix size = 256 x 256, 
echo time = 3.76 ms, repetition time = 7.99 ms, field of view = 
256 mm, slice thickness = 1.0 mm) acquired from a 3T Philips 
Gyroscan NT scanner (Royal Philips Electronics, Amsterdam, The 
Netherlands). The dataset included 13 healthy subjects (mean 
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age 10.42 years old; 5 males, and 8 females), 6 male subjects 
with autism spectrum disorder (ASD) (mean age 9.74 years old) 
and 11 subjects diagnosed with Attention Deficit/Hyperactivity 
Disorder (ADHD) (mean age 10.2 years old; 4 males, and 7 
females). The second dataset included sixteen whole brain, high 
resolution Tl -weighted 3D-volume MPRAGE images. Eight were 
healthy subjects (mean age 20.75 years old; all female) acquired 
from a Philips Intera 3T scanner (Royal Philips Electronics, 
Amsterdam, The Netherlands; matrix size = 256 x 256, echo 
time = 3.8 ms, repetition time = 8.1 ms, field of view = 256 mm, 
slice thickness = 1.0 mm). Four were diagnosed with primary 
progressive aphasia (mean age 68 years old) acquired from a 1.5T 
Philips Gyroscan NT (Royal Philips Electronics, Amsterdam, The 
Netherlands; matrix size = 256 x 256, echo time = 32 ms, repe- 
tition time = 6.85 ms, field of view = 230 mm, slice thickness = 
1.0 mm). The other 4 were subjects with Alzheimer's disease (AD) 
(mean 75.6 years old) acquired from a 1.5T Philips Gyroscan NT 
(Royal Philips Electronics, Amsterdam, The Netherlands; matrix 
size = 256 x 256, echo time = 3.2 ms, repetition time = 6.9 ms, 
field of view = 240 mm, slice thickness =1.2 mm). All data were 
resampled to isotropic 1 mm 3 voxel resolution and image size of 
181 x 217 x 181. 

LDDMM 

Given an atlas image Iq and a target image I\ , which can be rep- 
resented as functions Iq, C _R 3 — >■ R on the spatial domain 
Q C R 3 , the LDDMM algorithm (Beg et al., 2005; Ceritoglu 
et al., 2009) computes a diffeomorphic transformation cp:£2 — > 
Q. between these images such that I\ = Irj-cp . The diffeomor- 
phism cp = (f>i is defined as the end point of the curve cp — cp^, t € 
[0, 1] satisfying the ordinary differential equation cj> f = v f ((J) f ), 
where cpo = Id is the identity transformation and v t :Q. —> R 3 , t G 
[0, 1] is the time dependant velocity vector field of the flow of 
deformation. 

The diffeomorphism cp is calculated as: cp = cpi = v t (cp f ) dt 
with cpo = Id, where the optimal v t is estimated by solving the 
variational problem: 

v= argmin ( f \\Lv t \\ 2 2 dt + \\l 0 ■ cp^ 1 - h f \ (1) 

Theoretically, the minimizer of equation ( 1 ) results in a geodesic 
path for the curve cp t in the space of diffeomorphisms. To ensure 
that the solution lies in the space of diffeomorphisms, smoothness 
is achieved by defining the operator L as: L = — otV 2 + 7/3x3, 
where V 2 is the Laplacian operator. 

In LDDMM, steepest gradient descent approach is used to per- 
form the minimization in equation (1) and the velocity field at 
each gradient descent iteration k is updated with 

v? +1 = vf-s(V vf £ t ), (2) 

where V v £ t is the gradient of the cost in equation (1). 

V v £ f = 2v, - K * ( 4 | £<K, 1 1 V/ f ° (/ f ° - J*) Y (3) 



where cp s , t = cp f • cp s 1 ,/ t ° = I 0 • $t,oJ} = h ■ <k l, |-Dcp;' j is the 

determinant of the Jacobian matrix, K = (/J L) 1 and * is the 
convolution operation. 

In the numerical implementation of LDDMM, the time 
parameter t of the flow is discretized with a fixed total number 
of timesteps T, where T = 10 is selected as the default value in 
general. Selection of a smaller T causes gradient descent to termi- 
nate with a higher final mismatch error 1 7o • cp~ 1 — I\ | ^ 2 between 
the registered atlas image and the target image. In this article 
T = 1 is chosen to approximate the small deformation setup for 
comparison with T = 10. 

The background space of the images is represented with £2 = 
[0, l] 3 and the convolution operation in equation (3) is calculated 
in Fourier domain. The operator K acts as a low pass filter at each 
iteration of gradient descent and the parameters a and y controls 
the amount of smoothing and the elasticity of the deformation. 
Selection of these parameters depends on the size of the deforma- 
tion (in pixels) necessary to register the features of the atlas image 
to the features of the target image. Resampling the atlas and tar- 
get images to same resolution and doing an initial rigid or affine 
alignment as a preprocessing step for LDDMM allows the selec- 
tion of a/y in the range of 0.01-0.001 to accurately register the 
features of MR brain data with a typical resolution in the range 
of 0.5-2 mm/voxel. As a/y decreases in the 0.01-0.001 range, the 
image matching quality increases and the deformations become 
more local and elastic (Beg et al., 2005; Ceritoglu et al., 2009; 
TwardetaL, 2011). 

In the cascading implementation of LDDMM (Ceritoglu et al., 
2009), The final diffeomorphic deformation cp between images 
Jnand-fi is calculated by combining deformations cpi, cp2, . . . , cp„ 
and cp = cp^ 1 ■ cp^ 1 • ... • cp^ 1 where (i) cpiis calculated between 
Inand I\, (ii) cp2is calculated between Iq ■ cp^and I\, (hi) cp3 is 
calculated between Iq ■ cp^ 1 • cp^and I\. In each step, the regu- 
larization parameter y is fixed to be 1 and a is decreased to ensure 
that final registration is more robust and does not converge to an 
apparently incorrect solution compared to computing cp directly 
between Iq and I\ using very small a value (Ceritoglu et al., 2009) 
directly without cascading. The cascading implementation does 
not calculate a geodesic trajectory between the images and the 
mathematical properties of LDDMM such as metric distance and 
the encoding of the diffeomorphic flow of the shape with the 
initial tangent vector are not retained. 

It is worth noting that geodesic trajectories can be obtained 
via multi-scale approaches (Risser et al., 2011; Vialard et al., 
2012) where a weighted sum of Gaussian kernels is used to 
define the operator K in equation (3) control the elasticity of the 
deformation at each scale. 

SEGMENTATION VIA SINGLE-ATLAS 

The single-subject atlas is the JHU-DTI-MNI atlas (Oishi et al, 
2009), which is a single-subject template in the ICBM-DTI-81 
space. The template has an isotropic 1 mm 3 voxel resolution and 
image size of 181x217x181. Histogram matching was applied 
to the atlas and subject image, followed by initial affine align- 
ment of the atlas onto the subject via AIR (Woods et al., 1998a,b). 
MRIStudio was then used to register the subject to the atlas. The 
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resulting diffeomorphism was used to transfer the ROIs defined in 
the atlas to the subject space resulting in ROIs segmented in the 
subject space. This procedure is illustrated in Figure 1. All exper- 
iments were performed using Intel Xeon CPU E5530 at 2.40GHz 
with 16 cores. 

SEGMENTATION VIA MULTI-ATLAS FUSION 

The algorithm for multi-atlas segmentation using LDDMM (Tang 
et al., 2012, 2013) is briefly summarized as follows. Assume there 
are multiple MRI atlases, each of which contains a collection of 
locally-defined charts obtained from manually segmented struc- 
tures. A MAP approach is used to estimate the high dimensional 
segmentations from the class of generative models representing 
the observed MRI which is assumed to be a Gaussian random 
fields conditioned on the atlas charts and diffeomorphic change 
of coordinates of each chart. The charts and their diffeomorphic 
correspondences are unknown and viewed as latent or hidden 
variables. The expectation-maximization (EM) algorithm yields 
the likelihood-fusion equation which is maximized by the a pos- 
teriori estimator of the segmentation labels. The fused likelihoods 
are modeled as conditional Gaussian random fields with mean 
fields a function of each atlas chart under its diffeomorphic 
change of coordinates onto the target. The conditional-mean in 
the EM algorithm specifies the convex weights with which the 
chart-specific likelihoods are fused. The multiple atlases with the 
associated convex weights imply that the posterior distribution 
is a multi-modal representation of the measured MRI. As with 
the single-atlas segmentation, the resulting diffeomorphism was 
used to transfer the ROIs to the subject space resulting in ROIs 
segmented in the subject space. 

COMPARISON METRICS 

Automated segmentation of subcortical structures was compared 
with manual segmentations of the same subjects. The man- 
ual segmentations are based on anatomical definitions (Wang 
et al, 2007; Qiu et al., 2009b, 2010) with the assistance of an 
atlas (Mai et al., 1997); for details see http://caportal.cis.jhu.edu/ 
protocols. Manual segmentation was performed with the open 
source software Seg3D which follows the radiological convention 
for displaying images. 

Accuracy for each ROI was quantified by three metrics. 
The first is the kappa statistic (Cohen, 1960) is defined by 

K = (pagree ~ Prandom) / (l ~ Prandom) where p agree IS the fraction 



of voxels in which the given segmentation agrees with the 
manual segmentation, and prandom is the fraction one would 
expect by random chance (based only on the volumes of fore- 
ground and background). Note that ids biased by the volume 
size of the structure; generally, the bigger the structure, the 
higher the kappa statistic is. For brain structures, a value of 
k = 0.8is considered quite good. The second is the volume 
error which quantifies volume difference between two label 
defined by VD(L A ,L M ) = 100|V(L A ) - V(L M )\/V(L M ) where 
V (L A )and V (Lm) are respectively the volume size of the auto- 
mated and manual segmentation. The third is the L\ misclas- 
sification error (Miller et al, 2000) defined by L\ = \A U M| — 
| A n M| / | A U M|/2 where A and Mrespectively are the set of 
automated and manual labeled voxels. 

STATISTICS 

In the first experiment, One-Way ANOVA was performed on k 
and L\ from the automated and manual segmentations. Multiple 
comparisons were conducted using Tukey's HSD (honestly sig- 
nificant difference) test between pair wise groups. In the second 
experiment, two statistical analyses were performed to examine 
the effects of LDDMM parameters — T and a — on the segmenta- 
tion accuracy of single-atlas based registration. In the first analy- 
sis, meanii and volume errors were calculated for each subjects 
by averaging the errors of the segmented ROIs for that subject. 
Then Two-Way ANOVA was conducted to examine the effects of 
changing T and a on mean Iiand volume error. Multiple com- 
parison tests using Tukey's HSD was carried out after ANOVA. In 
the second analysis, the effect of T on the errors for each ROI was 
analyzed separately via one sided paired samples f-test between 
T = 1 and T = 10 results for given a. 

RESULTS 

For the first experiment, multi-atlas segmentation with cas- 
cading values of a = 0.01, 0.005, 0.002 and T = 10 was com- 
pared with single-atlas segmentation with the same values of 
a and T, FSL (version 5.0) and FreeSurfer (version 5.2.0). 
Six structures were analyzed (left and right pairs of cau- 
date, putamen and globus pallidus). Mean and standard devi- 
ation of kappa overlap and L\ errors are shown in Figure 2. 
One-Way ANOVA results (F L . caudate = 155.57, F R . caudate = 196.8, 

fL.globuspallidus = 95.45, FR.pallidus = 96.85, -FL.putamen = 182.73, 

FR.putamen = 228.23, df = 3, P< 0.001) for kappa overlap 

and (-F L . C audate = 163.08, FR.caudate = 187.79, FLglobuspallidus = 
73.89, -FR.globuspallidus = 75.90, FL.putamen = 186.88, FR.putamen = 

244.95, df = 3, P < 0.001) and the multiple comparisons test 
with Tukey's HSD showed that for each of six structures con- 
sidered, the automated segmentations obtained from multi-atlas 
LDDMM are statistically more accurate than those from the other 
three methods. For the second experiment, the following combi- 
nations were used: a = {0.01, 0.005, 0.002, 0.002with cascading} 
and T = {1, 10}. So, a total of (4 a-values) x (2 T- values) = 8 
different automated segmentations were calculated for each sub- 
ject using LDDMM single-atlas segmentation. Fourteen struc- 
tures (left and right pairs of caudate, putamen, globus pal- 
lidus, hippocampus, amygdala, thalamus, ventricle) were ana- 
lyzed. Figure 3 (top and middle rows) shows the mean and 




FIGURE 1 | Example axial slices of atlas and subject T1 weighted MR 
image and the boundaries of their ROI labels at each registration step. 

(A) Original atlas image and atlas ROIs; (B) atlas image and ROIs after 
linear AIR transformation; (C) atlas image and ROIs after AIR and LDDMM 
transformations; (D) subject image and subject ROIs. 
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LI Error 
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a = 0.005 a = 0.002 a =0.002 

(cascading) 



Volume Error 





I! Ill 



a = 0.005 a = 0.002 a =0.002 

(cascading) 



Run Time (sees) t 

J J J J 

a = 0.01 a = 0.005 a - 0.002 a = 0.002 

(cascading) 

FIGURE 3 | Mean and standard deviation of mean L q error (top), mean 
volume error (middle) and LDDMM running times (bottom) calculated 
from a population of 16 subjects for different values of a, T = 1 (blue) 
and T = 10 (red). Mean errors were calculated using 7 pairs of caudate, 
putamen, pallidus, amygdala, hippocampus, thalamus, and ventricle for 
each subject. 



Kappa Statistics 




L Caudate R. Caudate L. Pallidus R. Pallidus L Putamen P. Putamen 




L Caudate P.. Caudate L. Pallidus R. Pallidus L Putamen R. Putamen 



FIGURE 2 | Mean and standard deviation of mean Kappa statistics 
(top) and mean Li errors (bottom) for left and right pairs of caudate, 
putamen and pallidus in 30 subjects via four different segmentation 
methods: multi-atlas (blue), single-atlas (red), FSL (green), and 
FreeSurfer (purple). 



standard deviation of these errors computed over 16 subjects. 
There was no significant interaction between the effects of T 
and a, (F = 1.6785, df = 3, P = 0.17527). Simple main effects 
analysis showed that the effects of both a (F = 7.5801, df = 3, 
P = 0.00011) and T (F = 9.3323, df = 1, P = 0.0028) are sta- 
tistically significant. A similar analysis on mean volume error 
showed that there was no significant interaction between the 
effects of T and a, (F = 0.16, df = 3, P = 0.9230). Simple main 
effects analysis of the mean volume error showed that the effect 
of a (F = 4.1938, df = 3, P = 0.0073) was statistically significant 
but the effect of T (F = 2.7054, df = 1, P = 0.1026) was not sta- 
tistically significant. Multiple comparisons test with Tukey's HSD, 
showed that there is a significant difference (P < 0.01) between 
using a high a (0.01) and low a (0.002 with and without cas- 
cading). Tables 1, 2 shows the statistically significant differences 
{P < 0.01) for Li error and volume error respectively for the 
results at T = 1 and T = 10 at a significance alpha level of 0.01. 

Figure 3 (bottom row) shows the average running times of 
LDDMM calculations for any selection of a and T for the pop- 
ulaton of 16 subjects. Computational complexity of each gradient 
descent iteration in LDDMM depends on T linearly. Smaller a 
also increases the computation time. 

DISCUSSION 

In this paper, we examined the computational complexity in map- 
ping multiple atlases for segmentation of subcortical structures — 
caudate, putamen, and globus pallidus. Our results show that a 
cascade of a parameters and selecting fewer time steps can reduce 



computational complexity as much as five times while obtaining 
reliable segmentations. 

We also compared the segmentations with two other 
methods — FSL and FreeSurfer. LDDMM based segmentation 
is capable of achieving comparable or superior accuracy as 
measured by either the kappa statistic or the L\ error. Extending 
the single-atlas LDDMM based segmentation to multiple atlases, 
we found that the segmentation accuracy was significantly 
increased. The even split between male and female subjects 
in the dataset of 30 subjects suggests that our algorithm is 
robust with respect to gender. Then in terms of pathological 
variability, the datasets consisted of control subjects, subjects 
with autism, as well as subjects with ADHD. So our multi- 
atlas LDDMM segmentation algorithm can be applied to a wide 
variety of subjects. This is crucial for clinical neuroimaging stud- 
ies since the ultimate goal of any segmentation algorithm is 
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Table 1 | For each ROI and for each a in LDDMM, the L\ error results 
for 7=1 and T — 10 were compared using one sided paired sample 
f-test. 

Paired f-test results for L-\ error 

a - 0.01 a - 0.005 a = 0.002 a = 0.002 

(cascading) 



T1 > T10 T1 > T10 



T1 > T10 

T1 <T10 T1 >T10 T1 >T10 

T1>T10 T1>T10 T1>T10 T1 > T10 



Left amygdala 
Left caudate 
Left hippocampus 
Left globus pallidus 
Left putamen 
Left thalamus 
Left ventricle 
Right amygdala 
Right caudate 
Right hippocampus 
Right globus pallidus 
Right putamen 
Right thalamus 

Right ventricle T1 > T10 T1 > T10 T1 > T10 T1 > T10 

Number of ROIs with 3 4 8 3 

significant 
differences 

Statistically significant (p < 0.07^ error differences between large deformation 
map (T= 10) and small deformation map (T= 1! are indicated by either "T1 > 
T10" or "77 < no." 



T1 > T10 T1 > T10 

T1 > T10 
T1 > T10 



Table 2 | For each ROI and for each a in LDDMM, the volume error 
results for 7 = 1 and 7 = 10 were compared using one sided paired 
samples t-test. 

Paired t-test results for volume error 



a = 0.01 a = 0.005 a = 0.002 a = 0.002 

(cascading) 

Left amygdala 

Left caudate T1 < T10 T1 < T10 

Left hippocampus T1 < T10 T1 < T10 

Left globus pallidus 

Left putamen 

Left thalamus T1 < T10 T1 < T10 T1 > T10 

Left ventricle T1 > T10 T1 > T10 T1 > T10 T1 > T10 

Right amygdala 

Right caudate T1 < T10 T1 < T10 

Right hippocampus T1 < T10 T1 < T10 T1 < T10 

Right globus pallidus T1 < T10 T1 < T10 

Right putamen 

Right thalamus T1 < T10 

Right ventricle T1 > T10 T1 > T10 T1 > T10 T1 > T10 

Number of ROIs with 6 8 5 4 

significant 
differences 

Statistically significant (p < 0.01) error differences between large deformation 
map (T= 10) and small deformation map 77= 1) are indicated by either "T1 > 
T10" or "T1 < T10." 



to be able to accurately delineate brain structures with differ- 
ent diseases. The performance of multi-atlas LDDMM, in terms 
of segmenting subcortical structures of a different dataset, has 
been compared with STAPLE which is a classic label-fusion 
based segmentation method (Warfield et al., 2004). The seg- 
mentations were demonstrated to be superior to those from 
STAPLE (Tang et al., 2012). The high segmentation accuracy 
obtained from multi-atlas LDDMM is at the cost of more run- 
ning time, given that we need to first segment multiple atlases 
and we need to do multiple registrations. This issue led to 
the second experiment — analyzing the two parameters ot and 
T that determines the speed, memory requirement and accu- 
racy of LDDMM, trying to find the optimal a and T that 
will speed up the registration without affecting segmentation 
accuracy. 

In LDDMM, a and T are the two most important parameters 
that affect the smoothness of the computed diffeomorphism and 
the registration accuracy between the images that are registered. 
The results of the second experiment showed that decreasing a, 
decreased the mean and standard deviation of the errors and 
the results are statistically significant. Increasing T and using a 
large deformation setup instead of a small deformation setup 
also decreased the errors (Table 1) but results are less signifi- 
cant. Therefore, in general a small a should be selected for better 
registration. When each ROI is considered separately (first 3 
columns of Table 1 and Table 2), using a large deformation map 
(T = 10) instead of a small deformation map (T = 1) decreased 



the errors further with smaller a for given ROIs. Instead of using 
a single a to compute the final transformation with LDDMM, 
using a cascading scheme and decreasing a gradually to com- 
pute the final transformation allows the selection of a smaller T 
(last column of Table 2 and Table 2). This result is very impor- 
tant especially for multi-atlas LDDMM where a large number 
of LDDMM computations are necessary for superior segmen- 
tation accuracy. Thus in the extreme case of selecting T = 1 
leads to the mappings that capture only small deformations 
while increasing T with finer discretization of time flow leads 
to the mappings that capture larger deformations and are more 
smooth. 

Considered by some to be the "state of art" registration 
method (Ashburner and Friston, 2011), LDDMM has been 
shown to yield accurate segmentation of subcortical struc- 
tures in diverse datasets. Accuracy is obtained at a cost of 
computational complexity. As in the original pipeline (Miller 
and Qiu, 2009) in which segmentations were used to drive 
geodesic representations from which momentum and statis- 
tics were performed, the cascade approach described here 
allows us to generate intermediate analyses such as segmen- 
tations in a computationally efficient way. This is extremely 
important in our multi-atlas formulation with complex- 
ity 0(N), where N is the number of atlases. Thus via 
MRIStudio, LDDMM offers the potential for more precise 
and sensitive analysis of anatomical structures in neuroimaging 
studies. 
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